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Abstract 

We present a new numerical scheme for solving the advection equation and its 
application to Vlasov simulations. The scheme treats not only point values 
of a profile but also its zeroth to second order piecewise moments as depen- 
dent variables, for better conservation of the information entropy. We have 
developed one- and two-dimensional schemes and show that they provide 
quite accurate solutions within reasonable usage of computational resources 
compared to other existing schemes. The two-dimensional scheme can accu- 
rately solve the solid body rotation problem of a gaussian profile for more 
than hundred rotation periods with little numerical diffusion. This is cru- 
cially important for Vlasov simulations of magnetized plasmas. Applications 
of the one- and two-dimensional schemes to electrostatic and electromagnetic 
Vlasov simulations are presented with some benchmark tests. 

Keywords: Advection equation, Conservative form. Multi-moment, 
Information entropy, Vlasov simulations. Magnetized plasmas 



1. Introduction 

The kinematics of collisionless plasmas has been studied in a wide vari- 
ety of fields, such as in laboratory plasma physics, space physics, and astro- 
physics. Evolution of collisionless plasmas and self-consistent electromagnetic 
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fields is fully described by the Vlasov-Maxwell (or Vlasov-Poisson) equations. 
Thanks to recent development in computational technology, self-consistent 
numerical simulations of collisionless plasmas have been successfully per- 
formed from the first-principle Vlasov-Maxwell system of equations. 

There are two numerical methods to solve the Vlas ov equation. The most 
popu lar one is the Particle-In-Cell (PIC) method (IBirdsall and Langdonl . 



199ll ). which approximates the plasma by a finite number of macro-particles. 



Their trajectories calculated from the equation of motion are continuous in 
space, whereas electromagnetic fields are calculated on grid points in space. 
The PIC method has been used for a wide variety of plasma phenomena, 
because it gives satisfying results even with a relatively small number of par- 
ticles. However, the PIC method inherently has the large statistical noise 
due to an approximation of the distribution function by a finite number of 
particles. This noise only decreases in 1/ -y/iV when the number of particles 
N is increased, making it difficult to study such as particle acceleration and 
thermal transport processes, in which a small number of high energy particles 
play an important role. 

To overcome this problem, an alternative method free from the statistical 
noise has been used, in which the Vlasov equation is directly discretized on 
grid points in phase space. The so-called Vlasov simulation involves solving 
the advection equation in multidimensions (up to six). However, it has been 
widely known that a numerical solution of the advection equation suffers from 
spurious oscillations and numerical diffusion. A highly accurate scheme is 
required to preserve characteristics of the Vlasov equation (i.e., the Liouville 
theorem) as much as possible. In contrast to the PIC simulation, no standard 
scheme for the Vlasov simulation has been established thus far. 



Cheng and KnorrI ( 119761 ) proposed a splitting scheme, in which the elec- 



trostatic Vlasov equation is split into two advection equations in one-dimensional 
configuration and velocity spaces, and then are alternately advanced. Both 
equations reduce to a simple form of the linear advection equation. Moreover, 
their splitting method is equivalent to the second-order symplectic integra- 
tion so that the conservation of energy is very well. 

Following them, many authors have proposed high order advection schemes 
and applied them to the Vlasov simulation. One of standard time integra- 
tion methods for the advection equation is the semi-Lagrangian method, 
which advances a physical variable by interpolating its profile between grid 
poi nts and then following the c harac teristics backwar d in ti me. For exam- 
ple, ISonnendriicker et al.l (Il999[ ) and ICrouseilles et al.l (|2009[ ) employed this 
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method with the cubic B-sphne interpolation. iFilbet et al.l (l200l[ ) developed 
a semi-Lagrangian scheme in a conservative form with an upwind-biased La- 
grange polynomial interpolation, called the Positive and Flux Conserv ative 
(PFC) scheme. The PFC preserves mass and positivity. lUmedal (120081 ) pro- 
posed a non-oscillatory type of the PFC scheme. 

The above and many popular schemes consider the time integration of a 
single physical variable. On the other hand, the concept of " multi- moment" 



that t reats multiple dependen t variables has been proposed. lYabe and Aoki 



( I1991I ) and lYabe et al.l ( 1l99ll ) developed a multi-moment semi-Lagrangian 
scheme, called the Constrained Interpolation Profile (CIP) scheme. The 
CIP scheme employs the cubic Hermite interpolation. Difference of the CIP 
scheme from the others is that it treats not only a profile but also its first 
derivative in space as dependent variables, and their governing equations 
are solved as coupled equations. The CIP scheme provides comparable or 
better solutions even with a relatively smaller number of grid points com- 
pared to the others, although it requires a higher memory cost to store 
multiple dependent variables. Due to its high capab ility, the CIP scheme 
has b een applied not only to the Vlasov simulation fiNakamura and Yabel. 



1999) but also to rnagnetohydrodynamic si mulations ( Kudoh et al. . 19981 : 



Matsumoto and Sekil . I2OO8I : lYagi et al.l . l2009f l. 



Conservative schemes are preferable for the Vlasov simulation, because 
the integration of t he distribution fun ction in phase space is equ al to the 
number of particles. Yabe et al. ( 2001 ) and Takizawa et al.l (2002) also pro- 
posed a conservative form of the CIP, called the CIP-CSL2 (CIP-Conservative 
Semi-Lagrangian scheme with a second-order polynomial). The CIP-CSL2 
scheme treats point values of a profile and its cell-integrated values as depen- 
dent variables. The cell-integrated value is advanced in a conservative form to 
guarantee the conservation of mass. Therefore, the CIP-CSL2 scheme will be 
more suitable than the CIP scheme for the Vlasov simulation. Various CIP- 



CSL type schemes have been proposed thus far (e.g.. iTanaka et al.l. I2OOOI). A 



non-o scillatory type of the CIP-CSL scheme was proposed bv lXiao and Yabe 
(120011). and is su c cessfully extended t o fluid dynamic eq uations (IXiaol I2OO4J : 



Xiao et al.l . l2006t lli and Xiaol . l2007l ) . Ilmai et al.l (120081 ) developed a conser- 



vative Eule rian scheme, which is applied to the electrostatic Vlasov-Poisson 



simulation (jimadera et al.l . l2009l ) 



Although there has been a number of numerical schemes, they have been 
applied mainly to the electrostatic Vlasov-Poisson simulation. There is an 
increasing interest for the application of electromagnetic Vlasov simulations 
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to magnetized plasmas (IMangenev et al.l. l2002l: ISchmitz and Grauerl. l2006b : 
Sircombe and Arberl . l2009l : ISuzuki and Shigeyamal . I2OIOI : lUmeda et all . |2009| . 
2010l ). However, the Vlasov simulation of magnetized plasmas is more dif- 
ficult than the electrostatic one, because no suitable scheme for solving the 
gyr o motion (solid body rotatiori in velocity space) has been established thus 
far. ISchmitz and Grauerl fj2006al ) proposed a rigorous time-splitting method 
to solve the solid body rotation problem using the one-dimensional PFC 
scheme, called the backsubstitution. Even with this method, however, we 
find considerable numerical diffusion in a long time calculation that causes 
non-physical plasma heating during the gyration. For this reason, it is still 
quite limited to deal with plasma phenomena such as heating and accelera- 
tion on the basis of the electromagnetic Vlasov simulation. 

In this paper, we propose a new numerical scheme for the advection equa- 
tion, specifically designed to solve the Vlasov equation in magnetized plas- 
mas. We first argue in Section |2] that it is important to preserve high order 
moments of the distribution function to reduce numerical diffusion, on the 
basis of the concept of the conservation of the information entropy. We de- 
velop an advection scheme, in which not only point values of a profile but 
also its zeroth to second order piecewise moments are treated as dependent 
variables. Details of the scheme are described in Sections [3] and H] for one and 
two dimensions, respectively. Benchmark tests of the scheme and its appli- 
cation to electrostatic and electromagnetic Vlasov simulations are presented 
in Sections and El Finally, we summarize the paper in Section [71 



2. Concept of Multi-Moment Advection scheme 

We consider that the conservation of the information entropy is essential 
to develop a dissipationless scheme for the advection equation. When a 
profile f{x) follows the advection equation, Df/Dt = 0, its entropy function, 
h{x) = f{x)I{f{x)), also follows the equation. 



Dh f dl \Df 



where I{f{x)) is the self-information function of x, and f{x) is regarded 
as the probability function that rapidly decreases toward x = ±00 . The 
conventional information entropy function is h = —f log /. iTsallisI (jl988[ ) 
defined more general function of the form, h = {f — / (g — 1). The shape 
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of a profile will be preserved during the advection with constant velocity 
when its information entropy is conserved in a numerical simulation, 



h{x)dx = / f{x)I{f{x))dx = Constant, 

where Xi is the position in physical space. However, it is difficult to develop an 
advection scheme that exactly guarantees the entropy conservation, because 
I{f{x)) is nonlinear in general. Then, we approximate it by expanding in 
the Taylor series with respect to x, 



H 



E 



f 



/(/(0))+| 



x=0 



x" + ... 



x=0 



dx. 



This equation tells that the entropy can be described as a linear combination 
of the zeroth to n-th order moments, J {x'^f/n\)dx. We then consider that 
better conservation of the entropy may be achieved by preserving information 
of high order moments. For example, the entropy of a gaussian profile, f{x) = 
exp 



X — xq)'^ /2al]/ \/2Trax, is described by up to the second order moment. 



* = -E 
« E 

i 

= E 



f{x)\ogf{x)dx 



f{x) {x — Xo)^ dx 



fdx — 2xo 



xfdx + X fdx 

'Xi 



Since the gaussian profile is the most fundamental distribution in velocity 
space that maximizes the entropy in thermal equilibrium, preserving the 
entropy of this profile is crucially important for Vlasov simulations. We 
speculate from this discussion that it is essential to treat, at least, up to 
the second order moment as dependent variables, to keep the characteristics 
of the Vlasov equation. A "Multi- Moment Advection (MMA)" scheme is 
designed to treat the zeroth to second order moments as well as the profile 
itself as dependent variables, and advance them on the basis of their governing 
equations. 
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3. One-dimensional MMA scheme (MMAID) 

We consider the time evolution of a one- dimensional profile f{x,t) and 
its zeroth to second order moments defined as 

M"^ = —! x'^fdx, (m = 0,l,2). (1) 
m! J 

The conservative advection equation of / and governing equations of the 
moments are written as 

^ + u— = - — f (2) 

dt dx dx ' 

+ dx^ (uf) = 0, (3) 



dt J dx 

— / dx— {ux^f) = . — / ux'^-'fdx, (m = 1, 2) (4) 



dt ml J dx (m — 1)! 

where u is the velocity. Here, the conservative advection equation of / is cast 
into the advective form. Equations and are obtained by multiplying 
Equation ([2]) by x^/m! and then integrating over space. Note that the CIP- 
CSL2 scheme considers Equations ([2]) and ([3]) (m = only). Hereafter, we 
assume the constant velocity case, du/dx = 0, because we are concerned 
with the Vlasov equation. 

To solve a set of these equations, the one-dimensional MMA scheme 
(MMAID) treats four dependent variables; the point value of the profile 
fi, and the piecewise moments. 



1 



ml 



Mr+i/2 = —, / (m = 0,l,2), (5) 



and constructs a piecewise interpolation for / in a cell with a fourth order 
polynomial, 

5 

= ^kCk-i{x - Xi)^~^ , (6) 

k=l 

which gives an interpolation function for M™ as 

GTix) = 1 x"'FAx)dx 

5 

= ^A^(a;,a;OCfc;.(a;-a;0\ (m = 0,l,2), (7) 



k=l 
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where 



1, 

{kx + Xi) I {k 
| fc(fc+i) ^2 



A;XiX + xf}/{(A; + l)(/c + 2)} 



(8) 



To determine the coefficients C^-i, we use the variables at the upwind 
position as constraints, 

fit 

fiupj (9) 

= sgn(C)M-,,, 




[m 



0,1,2), 



where 



tup 
icell 



(10) 



i + sgn (C) , 
t + sgn (C) /2, 

are the position of the upwind grid and cell, sgn((^) stands for the sign of (, 
and ( = —uAt is the distance of the upwind departure position relative to 
Xi. Consequently, the coefficients are explicitly determined as 

fit 



Cl-i 
C2-i 



(11) 



4/i ~l~ fiup 



15sgn (C) 



2 

iup 



+ ixiXiup + x^i) Ml 



X 



-2 
Ax 
{{2x 

-2 {Axiup + 3x,) Ml^a + UMl^ii}] 





icell 



(12) 



C: 



3:i 
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2/. + /, 



lup 



2sgn(C) 
Aa;3 



X 



{(l0a;Lp + 25a;,a;„p + 7x2)M0,,, 
3sgn (C) 



(13) 



4;i 



Ax3 



4/^ + 3/, 



tup 



Ax^ 



■X 



{ (I5x2„p + 42xiXi,p + 13x2) Ml^ii 
-4 (18x,„, + 17x,) Ml^,, + 140M^,,}] 



(14) 



Ax4 



+ 



12sgn(C) 



X 



{( 



X, 



IMP 



5 [Xi. 



tup 



Ax3 

3xjXiup + ) M^^^ii 

+ x,)Mi,, + 10MLJ] 



(15) 
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where Ax = Xiup — Xj. 

Then, let us consider the time integration of the variables. The MMAID 
employs the semi-Lagrangian method. Equations (EJ-Q are split into two 
phases, advection and non-advection phases. We first calculate the advection 
phase as 

= "i^.(^. + C), (16) 
= "MlVi/2 + {"G°+i(.T,+i+C)-"G°(a:, + C)}, (17) 



a:i + l+C r^i+C, 



= "M™i/2 + {"Gr;i(.T,;+i +0 - "Gr(x. + C)}, (m = 1,2), (18) 

where the left-superscript indicates the number of time steps and the asterisk 
means that the variables are at the intermediate step. Here, we consider the 
case of the CFL number |nAt/Ax| < 1 for simplicity. We next advance the 
non-advection phase. Since u is constant, the right-hand side of Equation 
dl]) is equal to -uM™"^. Therefore, the non-advection phase is calculated as 

"+^^1/2 = *M^,/2 + n"+^MO.V2At, (19) 
"+iMti/2 = *M^V2 + ^(*^m/2 + ^"^'M°_i/2^) At (20) 

It is noteworthy that in the advection problem with constant velocity, the 
MMA scheme exactly guarantees the conservation of the zeroth to second 
order central moments defined as 

— / (x-ut)"'fdx, (m = 0,l,2). (21) 
ml 

The conservation of the zeroth order central moment (i.e., mass) is evidently 
guaranteed by Equation (fT7|) . Using Equations (fT8l) through (!20|) . the con- 
servation of the first order central moment is obtained as 

/oo 
{x - uAt) ■ ^+^fdx 
-oo 

i i 

POO 

= E"^m/2= / x--fdx, (22) 



and that of the second order central moment (i.e., variance) is also as 



1 /"°° 

-J (x - uAtf ■ "+^fdx 



i 

1 />oo 

= E"^m/2 = ^ / x^-Vtix. (23) 

4. Two-dimensional MMA scheme (MMA2D) 

We consider the time evolution of a two-dimensional profile /(x, y, t) and 
its zeroth to second order moments in the x and y directions defined as 

= JJ fdxdy (= M° = Mj) , (24) 
= ^ffx"^ fdxdy, (m = l,2), (25) 



m! 

= // (^ = 1'2). (26) 

The conservative advection equation of / and governing equations of the 
moments are written as 

+u^^ +v^^ — f^'^j^^'^\j ^27) 
dt dx dy \dx dy J ' 



dt J dx \m\ J J J dy \m\ 
= j j ux"'-' fdxdy, (m = 1,2), (29) 

+ / dx^ [ uy'^fdy) + [dy^(y^[ vfdx 



dt J dx \m\ J J J dy \m\ , 

. 1 1 vy"'-^ fdxdy, (m = 1, 2) , (30) 



where u and v are the velocity component in the x and y directions. Here, 
the conservative advection equation of / is cast into the advective form. 
Equations fl28|) - fl30|l are obtained by multiplying Equation fl27|) by x'^/m\ 
or y™'/m\ and then integrating over space. Hereafter, we assume du/dx = 
dv/dy = 0, but not necessarily du/dy = dv/dx = 0, and use vector forms 
= {M^, M^), X = {x, y), and u = {u, v). 
To solve a set of these equations, the two-dimensional MMA scheme 
(MMA2D) treats six dependent variables; the point value of the profile fij, 
and the piecewise moments, 

MT+i/2,j+i/2 = — / x"'fdxdy, (m = 0,l,2), (31) 

and constructs a piecewise interpolation for / in a cell with a polynomial, 

3 3 

Fidi^^y) = ^^lkCik;i,j{x - Xi)''"^ {y -yjf'^ , (32) 

1=1 k=l 

which gives an interpolation function for TVf" as 



=1 k=l 

where G"" = (G^, G"^), = = G^, and has been already defined in 
Equation (|8]). Note that the order of the interpolation function of / is same a s 



that used in the two-dimensional CIP-CSL2 scheme (ITakizawa et al.l . |2002| ). 

To determine the coefficients Cik-ij, we use the variables at the upwind 
position as constraints, 

Fij(^Xiiip, yj^ fiupjy 

'^iji.'^i^ Vjup) fijupy (2^) 

F^iji.-^iupy Vjup) fiupjupy 

G1^j{xiup,yj^p) = sgn (Oj) sgn (r/ij) (m = 0,1,2), 



where 



iup = i + sgn (Cij) , 

icell = i + sgn (Cij) /2, 

jup = j + sgn {r]ij) , 

jcell = j + sgn(?7ij) /2, 



(35) 
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are the position of the upwind grid and cell in the x and y directions, and 
iCijyVij) is the distance of the upwind departure position relative to {xi,yj), 
determined with a second order accuracy, 



-UjAt + Vi 
-ViAt + Uj 















dx y 





(Note again that we have assumed du/dx = dv/dy = 0, but not necessarily 
du/dy = dv/dx = 0.) Consequently, the coefficients are explicitly deter- 



mined, which are listed in Appendix A 



Then, let us consider the time integration of the variables. For the Vlasov 
equation in magnetized plasmas, we consider two problems of the advection 
with constant velocity and the solid body rotation. 

4.1. Time integration for the advection 

When the velocity is constant in space, it is easy to implement the semi- 
Lagrangian method to the two-dimensional scheme, because the volume ele- 
ment of the fluid does not change its shape. The advection phase of Equations 
(I271)-(l30l) is calculated as 



n+l r 



n+l j^jO 



icell,jcell 



"F,j(a;j + C,yj +77), 
sgn (C) sgn (r?) x 

Fij{x,y)dxdy 



(36) 



"M 



icelLjcell 



yj+TfJxi+C 

Vjup rXiup+C 

+ / "Fi^jup{x,y)dxdy 

ryjup+n f-xi^p+Q 
+ / "'Fiup,]up{x,y)dxdy 

J yjup J Xi^p 

sgn (C) sgn {rj) x 

/ x'^ ■''F,^j{x,y)dxdy 

yj+r]Jxi+C 
yjup i-Xiup+C 

/ x"" ■"F,up,j{x,y)dxdy 

yj+nlxi^p 



(37) 
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a;™ ■ "Fi^jupix,y)dxdy 



Xi+( 



+ 



Vjup+n rXi^p+C 



Vjup 



X 



■lup.jup 



{x, y)dxdy 



(m = 1,2). (38) 



Here, we consider the case of the CFL number < 1. Using Equation fl33|) . 
integrations in the right-hand side of Equations fl37|) and f l38|) are exactly 
calculated as, e.g.. 



Vr 

Vj+riJxi+C 



"Fi^j {x, y)dxdy 



Vj 



Vj+V fXi. 
Vj 



Vjup fXi+C 



Vj+V fXi+C. 



Vi 



''Fi^j{x,y)dxdy 



icelljjcell 



sgn (C) sgn (ri) 



We next advance the non-advection phase. Since it is constant, the right- 
hand sides of Equations (l29l) and (l30l) are equal to uM^~^ and vM^~^, 
respectively. Therefore, the non-advection phase is calculated as 



n+l 



n+1 



M 
M 



1 

j+l/2j+l/2 
2 

j+l/2j+l/2 



'^i+l/2j+l/2 + « 



At, 



-u 



i+l/2j+l/2 



-^i+l/2,j+l/2 + 



i+l/2j+l/2 



^^+1 /IfO 



(39) 



— At. (40) 



Similar to the one- dimensional case, in the two-dimensional advection 
problem with constant velocity, the MMA scheme exactly guarantees the 
conservation of the zeroth to second order central moments in the x and y 
directions. 

4-2. Time integration for the solid body rotation 

When the velocity varies in space, it is quite complicated to implement the 
semi-Lagrangian method, because the exact integration of the interpolation 
function over a distorted volume element is required. In such a general case 
(including the solid body rotation problem), it is rather convenier it to use the 
finite volume method for the time integration of the moments. Ili and Xiao 
(120071 ) implemented a convenient time integration method to the CIP-CSL2 
scheme, in which cell-integrated values are advanced by the finite volume 
method with the Runge-Kutta time integration, whereas point values are 
advanced by the semi-Lagrangian method. 
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The MMA2D employs their strategy. When the velocity is given as 
{u,v) = {—uy,ux), Equations (I251) - (15P1) are written in the following finite 
volume formulation, 



"™i+l/2,j + l/2 

Lodt 



(41) 



ujdt 



^^•^y"+l/2,j + l/2 

ujdt 



( [ yfdy) -( [ yfdy) 

-Iff xfdx) -( [ xfdx) \ , 

{\J /i+l/2,j \J /i+l/2,j + lj 

-±Lr(fyfdy) -xYlJfyfdy) 

— J ( [ x'^+^fdx^ - ( [ x^'+^fdx^ 

™! [\J /i+i/2,j \J Ji+l/2,j+l 

-^—(ffyx^^-'fdxdy) , (m = l,2), (42) 

(m- 1)! \J J /i+i/2,i+i/2 

~l( fy-^+^fdy] ~(fy"^+'fdy) 



/ xfdx 

1/2, j \J / i+l/2,j + l 



—L-^Q'J'xy^^-'fdxdy] , (m = l,2). (43) 



i+l/2,j+l/2 



Evidently, Equation (14T|) guarantees the conservation of mass. Line-integrated 
variables of / (up to third order) appearing in the right-hand side of Equa- 
tions (l4T!) - (H5l) are constructed from the interpolation function as 



— fx'^fdx) = sgn(Cij)^ /" x"^F,^j{x,yj)dx 

-I J icelLj ™- -Ixi 

3 

= sgn(C.J)^A™(a:™p,a:OC^lfc;^.,Ax^ (m = 0, 1, 2, 3) (44) 
fc=i 



sgn (77,.,) J2 ^riVjup, yj)Cn:;,jAy\ (m = 0, 1, 2, 3) , (45) 



1=1 



where we expand the definition of for the third order moment, 
A^(x,Xi) - 



(k + l) {k + 2) {k + 3) 

^3 _^ -v^ ^' XiX^ + kx\x ^x\ \ . (46) 



/c(fc + l)(fc + 2) 3 , k{k^\) 2 , , 2 , 3 
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The non-advection terms (third terms on the right-hand side) of Equations 
and are equal to — and for m = 1. For m = 2, the cross 
moment appears, which is also constructed as 

xyfdxdy] = sgn (C^j) sgn (77,^ ) / / xyFij{x,y)dxdy 

/ icelljcell Jyj J Xi 

= sgn(C.,j)sgn(7;,;j) X 

3 3 

)Cik-r^,Ax''AyK (47) 

(=1 fe=i 

For stable calculation, the time integration of Equations (I41D-(H31) is im- 
plem e nted with the third-orde r TVD Runge-Kutta method (jShu and Osherl . 



19881 : iGottlieb and Shul . Il998l ). To construct the variables (Equations 



[5|1 . and (H7|) ) at each Runge-Kutta stage, we should calculate the coeffi- 
cients of the interpolation function at each stage, which depend on / as well 
as M"^. Therefore, intermediate values of / at each stage are necessary. 
To obtain these, we first solve the equations of the characteristics with the 
Runge-Kutta method. 



dCi,j 
dt 

dt 



-UJ {yj + 7]ij) , with = 0, 
u {xi + Cij) , with °?7jj = 0, 



to get intermediate values of the upwind departure position at each stage, 
{xi + ^Cij^Uj + ^Vi,])^ where the left-superscript k denotes the Runge-Kutta 
stage. The solution is exactly obtained as 

+ = COS (— *^aa;At) Xi — sin i^—^auj^t) yj, 
Uj + ''Vij = sin (^—''aujAt) Xi + cos (— ^awAt) yj, 

where is the Runge-Kutta coefficient, (^a,^a) = (1,0.5) for the third- 
order TVD method. Thus, the intermediate values of / are calculated with 
the semi-Lagrangian method, 

fi,j Fiji^Xi + CiJiVj ~l~ Vi,j)- 

Consequently, we can calculate the coefficients at each stage, and then ad- 
vance the moments with the third-order TVD Runge-Kutta method, which 
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is implemented as follows, 



n+i_^m _ i'^M'" + -{2M"* + i?(2/,^M'")a;At}, (m = 0, 1, 2)(48) 
3 3 

where R stands for the right-hand side of Equations pT|) - fj43|) . Finally, we 
calculate the point value as 

^~^^fi,j = "-Pij(cos {-LoAt) Xi - sin {-ujAt) yj,sm {-ujAt) Xi + cos {-ujAt) yj).{49) 



Unlike the advection problem, the first and second order central moments 
are not necessarily conserved quantities in the solid body rotation problem. 
On the other hand, the sum of the square of the first order moments (i.e., 
the radius of rotation). 



\ 2 / ^oo POO \ 2 

xfdxdy ) + ( / / yfdxdy 

-ooJ —oo / \J~ooJ—co 

and the sum of the second order moments, 

r"00 poo 

(x^ + y^) fdxdy, 



(50) 



(51) 



oo«/ — oo 



should be conserved. Equations f H2]) and f H5]) (m = 2) guarantee the conser- 
vation of Equation fl^ . Taking the sum of Equations fH^ and fH51) (m = 1) 
over and then implementing the Runge-Kutta calculation, we obtain 



1 - 



""Ml - 00 At 



''Ml + ujAt 



{ujAtf 



(uAty 
6~~ 



Mi 



This is a third-order accurate solution of the centroid position of a profile. 
Equation flSUj) is numerically calculated as 

2 / \ 2 



15 



12 



36 



(52) 



Therefore, the MMA scheme will preserve the orbit of rotation with high 
accuracy, although it does not guarantee the exact conservation. 



5. Numerical tests of MMAID 

We perform some numerical tests with the MMAID, and compare them 
to the CIP-CSL2 and PFC schemes. The accuracy o f the C IP-CSL2 scheme 



is almost same as that of the GIF scheme flYabe et al.l . 120011). Comparisons of 



the CI P and PFC schemes with others are found in, e.g.. lFilbet and Sonnendriicker 



(2003 


), and 


Umeda 


( 


2008) 



to accurately give the initial condition of the piecewise moments. In the fol- 
lowing simulations, we determine the initial condition of the moments by a 
Gauss quadrature method. 

5.1. Advection problem 

We test the one-dimensional advection problem with constant velocity. 



dt dx 



of a gaussian profile. 



f{x,t = 0) = exp 



X 



2al 



(53) 



(54) 



with u = 1 and = 0.01. Since the numbers of dependent variables are 
different among the three schemes (four for the MMA, two for the CIP-CSL2, 
and one for the PFC), we use different grid sizes (0.01 for the MMA, 0.005 
for the CIP-CSL2, and 0.0025 for the PFC) so that the total memory usage 
is equal. The CFL number is 0.2. Figure [T] shows the results at t = 2.0 
solved by the MMA, CIP-CSL2, and PFC schemes. While the CIP-CSL2 
and PFC schemes show numerical diffusion, the MMA scheme provides quite 
an accurate solution. 

From simulation runs with different grid sizes, we evaluate the order of 
accuracy of the schemes. Figure [2] shows the local error at the center of the 
profile, 1/(0,0) — /(ut,t)|, as a function of the grid size. The slope of a line 
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corresponds to the order of accuracy in space. While the CIP-CSL2 (dia- 
monds) and PFC (crosses) schemes show the third order accuracy (dashed 
hne), the MMA scheme (triangles) achieves the seventh order accuracy (dot- 
dashed line). In Figure [2l the simulation results with the same memory usage 
are denoted as symbols with the same colors. The error of the MMA scheme 
is much smaller than the others with the same memory usage. 

Figure [3] shows the numerical test of the square wave advection after 4,000 
time steps. The grid size and CFL number are 0.005 and 0.2. Compared to 
the CIP-CSL2 and PFC schemes, the MMA scheme accurately keeps the 
discontinuity, although it shows small numerical oscillations. 

5.2. Electrostatic Vlasov simulations 

We apply the MMAID to electrostatic Vlasov-Poisson simulations. The 
one-dimensional electrostatic Vlasov-Poisson system of equations is written 
as follows, 

^ + ^,^_^^-0 (55) 
dt dx m dv ' 

^ ^ Aimoe ( [ fdv-l] , E = (56) 



where f{v,x) is the electron phase space distribution function, v is the elec- 
tron velocity, x is the position in the configuration space, e is the charge, 
m is the electron rest mass, rio is the density of ambient immobile ions, and 
E{x) and 0(x) are the electrostatic field and potential, respectively. 

To solve these equations, we treat eight dependent variables; point values 
of the distribution function, piecewise moments in the velocity space, and 
their cell-integrated values in the configuration space, 

fi,j f i^iy 

^mA, = ^ r\rf{v,x,)dv, 



ml 



/ij+1/2 = / f{vi,x)dx, 



= — / / V'^f{v,x)dvdx, 



ml 



where the subscripts i and j denote the grid position in the velocity and 
configuration spaces, and m = 0,1,2. We then split the Vlasov equation 
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(133]) into two phases, the advections in the velocity and configuration spaces, 
which are solved by the MMA and CIP-CSL2 schemes as follows: 

(1) The advection equation in the velocity space and its integration over 
the configuration space are 



dt m dv 
dfj+i/2 eEj+i/2 dfj+i/2 



0, 

0, 



dt m dv 

where we approximate j^''.^^ Efdx ^ -£'j+i/2/7+i/2- These equations are 

solved by the MMA scheme for {f,M"^) and {f,M"^), respectively. 
(2) The advection equations in the configuration space, 

dt dx ^ dt dx 

are solved by the CIP-CSL2 scheme for {f,f)- Multiplying these equations 
by v^/m\ and then integrating over velocity gives the governing equations 
of the moments. 

Since the (m+l)-th order moment is required for the advance of the m-th 
order one, the system is not closed (the closure problem). In this study, we 
approximate the equations as 

aM-,/, 9M-,/2 , ^ 

7^7 + 7^ = 0, — + Vi+i/2 7^ = 0. 

dt dx dt dx 

These advection equations are solved by the CIP-CSL2 scheme for (M°, M°), 
{M\M^), and {M^,M^). 

The Poisson equation fl56p is solved by the Fourier transform. The elec- 
tric charge density on the right-hand side can be immediately calculated by 
taking the sum of M° or M^/Ax over the velocity space. Since the MMA 
scheme guarantees the conservation of M° and M° (see Equation f lT7|) ). and 
the CIP-CSL2 scheme guarantees the conservation of / and M"^, our sim- 
ulation exactly conserves M°. Therefore, we use for the charge density 
calculation. 
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Physical variables in the system are normalized as follows; velocity by the 
thermal velocity Vth-, time by the inverse electron plasma frequency = 
■\/ m / A-Kn^^e^ , and position by the Debye length \d = Vth/^p- The boundary 
conditions are periodic in the configuration space, and open in the velocity 
space where constant incoming fluxes are assumed while outgoing flu x es are 



perfectly lost. We employ the splitting method of ICheng and KnorrI (jl976[ l 
for the time advance of the system. 

5.2.1. Linear Landau damping 

We first test the linear Landau damping. The initial condition of the 
electron distribution function is given as follows, 

/(u, X, t = 0) = ^= exp f-— ) [1 + 0.01 cos (0.5x)] . (57) 



27r V 2 

The simulation domain is [—4.5, 4.5] in the velocity space with 32 grid points, 
and [0,47r] in the configuration space with 32 grid points. The time step is 
0.01. FigureHlis the time profile of the electric field (normalized by the initial 
value), showing the exponential decrease. The simulation (solid line) agrees 
well with the linear theory (dashed line). The recurrence effect occurs at 
t ~ 42, which is comparable to the theoretical value t = 44.68 obtained from 



the free-streaming equation (jCheng and KnorrI . Il976l ) 



5.2.2. Two stream instability 

We next test the two stream instability. The initial condition of the 
electron distribution function is given as follows, 

f{v, X, t = 0) = -^v"^ exp (-^\ [1 + 0.05 cos (0.5x)] . (58) 

The simulation domain is [—5.0, 5.0] in the velocity space with 32, 64, or 128 
grid points, and [0,47r] in the configuration space with 64 grid points. The 
time step is 0.01. For comparison, we also perform the simulation in which 
the advections in the velocity space as well as the configuration space are 
solved by the C1P-CSL2 scheme. 

Figure |5] shows the electron distribution function at t = 25 and t = 50 
simulated with the MMA and C1P-CSL2 schemes. Both simulations show 
the well-known electron hole structure in the phase space. Figure |6] shows the 
velocity space distribution taken at x = 4.0 and t = 50. The MMA scheme 
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(blue line) describes finer structures than the CIP-CSL2 scheme with the 
same grid points (dashed line), and shows a similar result to the CIP-CSL2 
scheme with the double number of grid points (solid line). 

In Figure [71 we check the conservation of the total energy and the entropy, 



H = jj f{l-f)dvdx, 



with the different schemes (MMA or CIP-CSL2) and the different numbers of 
grid points (32, 64, or 128) in the velocity space. The MMA scheme provides 
the energy and entropy conservations well, even with the smallest number of 
grid points (32, green lines). It is noteworthy that the MMA scheme with 
32 grid points shows a better result than the CIP-CSL2 scheme with 128 
grid points (solid lines), which uses the memory twice as much as the MMA. 
Figure [S] shows the relative error of the total energy {e(t) — e(0))/e(0) at 
t = 75 obtained from the both schemes, as a function of the number of grid 
points in the velocity space. With the sufficiently large number of grid points, 
the both schemes show the same overall accuracy, because the same method 
is applied to solve the advection in the configuration space. 

6. Numerical tests of MMA2D 

We perform some numerical tests with the MMA2D, and compare them 
to the CIP-CSL2 and PFC schemes. We also apply the MMA2D to electro- 
magnetic Vlasov-Maxwell simulations and compare them to the PIC method. 

6.1. Solid body rotation problem 

We test long time simulations of the two-dimensional solid body rotation 
problem, 

^_y^ + ^^=0 (59) 
dt dx dy 



of a gaussian profile, 

/ (x, J/, t = 0) = exp 



{x - Xof {y - yoY 



2al 2al 



(60) 
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We co mpare the simulation s with the MMA to t he two-dimensional C IP 



CSL2 (ITakizawa et al.l . 120021 ) and backsubstitution ( ISchmitz and Graueii . l2006af ) 
schemes. Since the numbers of dependent variables are different among these 
three schemes (six for the MMA, four for the CIP-CSL2, and one for the back- 
substitution), we use the different numbers of grid points (34 x 34 for the 
MMA, 42 X 42 for the CIP-CSL2, and 84 x 84 for the backsubstitution) so 
that the total memory usage is equal. The simulation domain is [-0.5,0.5] 
in both the x and y directions. The open boundary condition is employed 
where constant incoming fluxes are assumed while outgoing fluxes are per- 
fectly lost. The time steps are 0.0047r for the MMA and CIP-CSL2 schemes, 
and 0.0027r for the backsubstitution scheme so that the CFL number is close 
among the three simulations. 

Figure |9] shows the results for a symmetric gaussian profile with xo = 
yo = 0, ax = (Jy = 0.1 after 50 and 300 rotations. The CIP-CSL2 scheme 
shows the most serious numerical diffusion. The backsubstitution scheme is 
better than the CIP-CSL2, but also shows considerable diffusion after about 
a hundred of rotations. On the other hand, the MMA scheme provides quite 
an accurate solution with little diffusion even after hundreds of rotations. 

Figure [10] shows the results for an asymmetric gaussian profile with xq = 
yo = 0,ax = 0.1, ay = 0.15. The MMA scheme provides a better solution 
than the others. The gaussian centroid (xo,|/o) and the standard deviation 
((Ta;, ay) calculated by the MMA scheme are kept correct even after hundreds 
of rotations, although a small error is seen on the tail of the profile (the 
outermost contour). It is also found from this simulation that the phase 
error of the rotation is very small. 

We next test the solid body rotation with the advection problem, 

f -(y-l/o)|f + (x-Xo)|^ = 0, (61) 

which describes the rotation around {x,y) = {xo,yo). For magnetized plas- 
mas, this corresponds to the electric field [E x B) drift motion. To solve 
this equation, we split it into two phases, the rotation phase (Equation (159|) ). 
and the advection phase, 

df ^ df df 
at ox oy 

Figure [11] shows the results for a symmetric gaussian profile with xq = 
—0.05, yo = 0.1, a^ = ay = 0.1. In the problem of the rotation with advection. 
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the backsubstitution scheme shows the most serious diffusion after several 
tens of rotations. Again, the MMA scheme provides a better solution with 
keeping (xq, yo) and {ax, (Jy) constant. Figure [T2l shows the temporal variation 
of {(Tx,(Ty) obtained by fitting the profile with the gaussian function. This 
result indicates that the MMA scheme can accurately solve the E x B drift 
motion for a long time with little numerical heating. 

From these simulation runs with different grid sizes, we evaluate the or- 
der of accuracy of the schemes. Figure [13] shows the error \f{x,y,t) — 
/(x, y, 0)\/N after 50 rotations, as a function of the grid size (iV is the number 
of grid points). All schemes show the third order accuracy in space (dashed 
line). In Figure [T^ the simulation results with the same memory usage are 
denoted as symbols with the same colors. The error of the MMA scheme 
(triangles) is ~ 10~^ times smaller than the others with the same memory 
usage. The errors of the schemes increase with time. 

We note that the order of the interpolation function of the MMA scheme 
(Equation fl32|) ) is exactly same as that of the C1P-CSL2 scheme. Neverthe- 
less, significant differences are found between the two schemes, suggesting 
that selection of constraints for the interpolation function is important to 
improve the accuracy. In this respect, we consider that the MMA scheme 
has the reasonable selection on the basis of the concept of the entropy con- 
servation (Section [2]), as it is designed specifically for Vlasov simulations. 

6.2. Electromagnetic Vlasov simulations 

We apply the MMA2D to electromagnetic Vlasov-Maxwell simulations. 
The one- dimensional electromagnetic Vlasov-Maxwell system of equations is 
written as follows, 

0, (s=p,e), (62) 

(63) 
(64) 
(65) 

where E and B are the electric and magnetic fields, j is the current density, 
c is the speed of light, is the charge, and the subscript s denotes particle 



— + n— + — 

dt ^ dx m. 



E 



dE 

'dt 
OB 

~dt 



cV X B — 47rj, 



-cV X E, 



3 = ^<ls j vfsdv, 



V X B 



dv 



s=p,e 
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species (p for protons and e for electrons) . We assume the two dimensionality 
in the velocity space, v = {v^, Vy, 0), E = {E^, Ey, 0), and B = (0, 0, B^). 

To solve these equations, we treat twelve dependent variables for both 
protons and electrons; point values of the distribution function, piecewise 
moments in the velocity space, and their cell-integrated values in the config- 
uration space, 

fi,j,k f i'^x,iy '^y,j J -^k^ 

MT+i/2j+i/2,k = ^1 / v"'f{v^,Vy,Xk)dv^du 



yi 



y-j 



fi,j,k+l/2 = / f{v^,i,Vyj,x)dx, 
-^i+i/2.,-+i/2,fc+i/2 = — T / / / v"'f{v^,vy,x)dv^dvydx, 



where the subscripts i, j, and k denote the grid position in the v^, Vy, and 
X directions, M'" = (M-,M-), = = M", M"" = {MZ^M:^,), 

= = M°, and m = 0, 1, 2. We then split the Vlasov equation 
into three phases, the advection in the c onfiguration space, the acceleration 



sp 

by electric fields, and the gyro motion (iGhizzo et all l2003l : lUmeda et al. 



20091 . |2010| ). The acceleration and the gyro motion are solved by the MMA 
scheme (Sections 14. Il and l4.2l) for (/, M"^) and (/, M ). The advection in the 
configuration space is solved by the C1P-CSL2 scheme for (/, /), (M", M"), 

{M\M ), and {M^,M ). Those are carried out in the same manner as the 
electrostatic Vlasov simulation that avoids the closure problem (see Section 



We solve the Maxwell equation with the implicit scheme (IHoshind . 11986 



19871 ). The current density can be immediately calculated from the sum of 



iVf-*- or M /Ax over the velocity space. We use the latter for the current 
density calculation, because they are conserved in the configuration space 
(but not in the velocity space). 

Physical variables in the system are normalized as follows; velocity by 
the speed of light, time by the inverse electron plasma frequency = 
y^m^T^T^e^; electromagnetic fields by an ambient magnetic field, and po- 
sition by the Debye length or the electron inertia length d^ = c/upe 
depending on problems. The boundary conditions are periodic in the con- 
figuration space and open in the velocity space. The time advance of the 
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system is implemented as follows; (1) advection in the configuration space 
with a half time step At/2, (2) advance of the electromagnetic fields with a 
half time step, (3) acceleration by electric fields with a half time step, (4) 
gyro motion with a full time step At, (5) acceleration by electric fields with a 
half time step, (6) advection in the configuration space with a half time step, 
(7) current density calculation, (8) advance of the electromagnetic fields with 
a half time step. This algorithm is a second order leap-flog time integration, 
also employed in PIC simulations. 

6.2.1. Perpendicular wave propagation 

We first test the wave propagation perpendicular to magnetic fields. The 
initial condition of the particle distribution function is given as follows. 



where Vg-th and ns{x) are the thermal velocity and the density distribution 
of s-species particles. We set a uniform density for both protons and elec- 
trons, and then add a small (1%) uniform random noise only for the electron 
density. The initial magnetic field distribution is uniform, Ey is zero, and 
Er^ is determined from the Poisson equation (!56|) . Simulation parameters 
are mp/nie = 16, Ve-th = 0.1, Vp-th = 0.025 (corresponding to a ratio of the 
proton to electron temperature being unity), and a ratio of the electron gyro 
to plasma frequency Uge/ojpe = 0.5. The simulation domain in the velocity 
space is [—4:Vs;th,4:Vs-th] with 50 grid points in both the Vx and Vy directions. 
The position is normalized by Ad. The grid size in the configuration space 
is equal to A^, and the spatial length is 512A/). The time step is 0.05/-\/2. 

Figure [T^ a) shows the Fourier spectrum of E^ obtained from the electro- 
magnetic Vlasov simulation integrated until t = 723.4. For comparison, we 
also perform the electromagnetic PIC simulation with the same parameters, 
which is shown in Figure [MT b). The number of particles is 5,000 in each cell 
so that the total memory usage is comparable between the two simulations. 
We can clearly identify the electron cyclotron (Bernstein) modes. We also 
identify the high-frequency X-mode and Z-mode at very low wavenumbers 
and frequencies close to the R- and L-mode cutoff (w = 2.56uge and 1.56uge), 
respectively. As expected, the Vlasov simulation provides a result with less 
noise, because it is free from the statistical noise. 

At very low wavenumber and frequency in Figure [T^ a). the contribution 
of protons is slightly seen. Then, we continue the same simulation for a long 



fs{Vx,Vy,X,t = 0) 




exp 




(s =p,e), 



(66) 
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time until t = 1447.4. The ion Bernstein modes and the lower-hybrid waves 
can be seen in Figure [T3 Even at this moment (electrons gyrate more than 
a hundred times), the total energy is conserved very well within an error of 
10~^, due to the fact that the MMA scheme can solve the solid body rotation 
problem with little numerical diffusion. 

6.2.2. Harris current sheet equilibrium 

We next test the one-dimensional Harris current sheet configuration, in 
which a dense plasma is confined at the center of the anti-parallel magnetic 
field configuration. The initial conditions of the particle distribution function 
and the magnetic field are given as a double Harris current sheet, for utilizing 
the periodic condition in the configuration space. 



fs{V^,Vy,X,t = 0) 
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(69) 
(,70) 



where uf is the diamagnetic drift velocity, Ugs = \qs\B()/msC is the gyro 
frequency outside the sheet, A is the half thickness of the current sheet, 
and L is the spatial length. The one-dimensional Harris current sheet is 
a steady-state solution of the Vlasov-Maxwell system, in which the plasma 
and magnetic pressures balance with each other. Therefore, it is a suitable 
benchmark test for electromagnetic Vlasov simulations: An accurate scheme 
is required for solving the advection and rotation in the presence of spatially- 
inhomogeneous plasma and magnetic field distributions and thereby keeping 
the equilibrium. 



Simulation parameters are mp/rrie = 10, Uge/ujpe = 0.1, v, 



e;th 



0.0707, 



Vp;th = 0.0224 (corresponding to a temperature ratio of unity), X = dp = 
c/upp (proton inertia length), and L = 1024A/). The thermal velocity is de- 
termined so as to satisfy the pressure balance. The position is normalized 
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by d^. The simulation domain is [—dvs-th, ^Vg-th] in the velocity space with 32 
grid points in both the and Vy directions, and [—L/2, L/2] in the configu- 
ration space with 1024 (Ax = Ad) or 256 (Ax = AXd) grid points. The time 
step is 0.027r. 

Figure [TBT a) shows the magnetic field and density (rie + Up) distribu- 
tions at t = 6283.19 (corresponding to 100 gyrations of electrons outside the 
sheet), obtained from the simulation with Ax = Xo- The simulation keeps 
the equilibrium of the current sheet with small numerical dissipation. The 
amplitude of numerically-produced electric fields is smaller than 10"'^. 

We check the pressure balance along the x direction between the magnetic 
field and plasma, 

where Pxx and Ts-^xx are the x-diagonal component of the pressure and temper- 
ature tensors. In the simulation with the MMA scheme, diagonal components 
of the pressure tensor can be easily estimated from dependent variables, 

= ^ [2M^^.._^i/2j+l/2 ~ 2V3,6Mj^.._^^/2,j+l/2 + ^xb^i+l/2,j+l/2] ) 
hj 

where Vxb is the bulk velocity in the x direction, also estimated from them, 

^ ^ // vJdVxdVy ^ E.J ^4;i+l/2,j + l/2 
JJ fdVxdVy J2ijK+l/2,j + l/2 ' 

Figure [TBT b) shows the pressure distribution. The pressure balance is pre- 
served very well within an error of 0.1 %. The plasma pressure is kept correct 
within an error of 2 %. 

For comparison, we also perform the simulation in which the CIP-CSL2 
scheme is applied to the advections in the velocity space as well as the con- 
figuration space. The number of grid points in the velocity space is 40 x 40 
for this simulation so that the total memory usage is comparable to the sim- 
ulation with the MMA scheme. The result is shown in Figure [T71 Compared 
to the MMA scheme, the CIP-CSL2 scheme causes the dissipation of the cur- 
rent sheet (increase of the sheet thickness), through the numerical heating of 
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the plasma. The proton and electron pressures increase by 10 % and 40 % 
around the edge of the sheet, respectively. The total pressure increases by 2 



70. 



Figure [18] shows the result from the same simulation with the MMA 
scheme, but with a larger grid size of Ax = 4Ad. The simulation again 
keeps the equilibrium of the current sheet. The order of a numerical error 
is comparable to the result with Ax = Xo- The result confirms that our 
electromagnetic Vlasov simulation code is stable even with the grid size larger 
than the Debye lengt h. This is an advantage ove r an explicit PIC simulation, 
as suggested by, e.g., lUmeda et al.l (120091 12010|). 



7. Summary and Discussion 

We have presented a new numerical scheme for solving the advection 
equation and the Vlasov equation. The present scheme solves not only point 
values of a profile but also its zeroth to second order piecewise moments 
as dependent variables, for better conservation of the information entropy. 
We have developed one- and two-dimensional schemes, and have shown their 
high capabilities. The present scheme provides quite accurate solutions even 
with smaller numbers of grid points, although it requires a higher memory 
cost than other existing schemes. We have shown that, however, the total 
memory usage of the present scheme is smaller than the others for the same 
accuracy of solutions 

Applications of the one- dimensional scheme to the electrostatic Vlasov 
simulations (linear Landau damping and two stream instability), and the two- 
dimensional scheme to the electromagnetic Vlasov simulations (perpendicular 
wave propagation and Harris current sheet equilibrium) have been presented. 
The two-dimensional scheme allows us to solve the gyro motion and the ExB 
drift motion for a long time with little numerical heating. Since the present 
scheme treats the zeroth to second order moments and advances them on 
the basis of their governing equations, the particle momentum and energy as 
well as mass are conserved very well. This is important for studying plasma 
phenomena such as convection, heating, and acceleration. 

Although the present scheme correctly solves the moments up to the sec- 
ond order, the entropy is numerically increased in the two stream instability 
(Section I5.2.2p . by the dissipation of fine structures in velocity space. This 
is understood that the perturbation of lower order moments decreases and 
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higher order moments contribute to the entropy w hen the so-called filamen 



tation phenomena proceed (e.g., iPohn et al.l . 120051 ). Since the filamentation 



is inevitable as long as discretizing velocity space, it is essentially impossible 
to exactly conserve the entropy with finite information. 

The present scheme is designed specifically to solve the advection and 
rotation in velocity space in the Vlasov equation. In our Vlasov simulation 
code, the CIP-CSL2 scheme has been employed to solve the advection in 
configuration space. We note that, however, another scheme can be applied 
to the advection in configuration space and be combined with the present 
scheme. 

Several advection schemes proposed for Vlaso v simulations are de s igned to 
prese rve positivity and non-oscillatory property (IFilbet et al.l . 1200 ll : lUmedal . 
20081 ). in order to suppress a non-physical growth of plasma waves caused 



by nu merically-produced positive gradient in velocity space. lUmeda et al 
( 120061 ) argued that the positivity-preserving and non-oscillatory are impor- 
tant properties for reliable Vlasov simulations. On the other hand, it may 
cause a non-physical global evolution driven by plasma pressure gradient 
that is increased by the numerical diffusion in velocity space, even when the 
plasma should be in an equilibrium state. Although the present scheme is not 
positivity-preserving or non-oscillatory, we have shown that the obtained re- 
sults are better than the others. We thus consider that preserving high order 
moments is another important property for Vlasov simulations. 

Magnetized plasma phenomena that allow to assume the two dimen- 
sionality in velocity s pace are limited, e.g., strictly perpendicular shocks 
( lAmano and Hoshind. 120091) and the two-dimensional Kelvin - Helmholtz in- 
stability (IMatsumoto and Hoshind . l2006l : lUmeda et al.l . |2010[ ). Even consid- 
ering the one dimension in configuration space, there are many phenomena 
that should treat the full three-dimensional velocity space. We are now 
developing a three-dimensional scheme and its application to the full electro- 
magnetic Vlasov simulation, which will enable us to study a wide variety of 
coUisionless plasma phenomena. 
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Appendix A. CoefRcients of the interpolation function of MMA2D 
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Figure 1: Numerical tests of the one-dimensional advection of a gaussian profile with (a) 
the MMA, (b) the CIP-CSL2, and (c) the PFC schemes. Solid lines with symbols show the 
simulation result, and dashed lines are the exact solution. In (a), the solid line overlaps 
the dashed line. The grid sizes are 0.01 for the MMA, 0.005 for the CIP-CSL2, and 0.0025 
for the PFC. 



35 




0.001 0.010 
Grid size 



Figure 2: Local errors of the one-dimensional adveetion of a gaussian profile as a function 
of the grid size. Triangles, diamonds, and crosses are the results obtained from the MMA, 
CIP-CSL2, and PFC schemes. Dashed and dot-dashed lines represent the third and sev- 
enth order accuracies. Symbols with the same colors are obtained from the simulations 
with the same memory usage. 
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Figure 3: Numerical tests of the one-dimensional advection of a square wave profile with 
(a) the MMA, (b) the C1P-CSL2, and (c) the PFC schemes. Solid lines with symbols show 
the simulation result, and dashed lines are the exact solution. 
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Figure 4: Time profile of the electric field (normalized by the initial value) in the linear 
Landau damping simulation with the MMA scheme. A solid line shows the simulation 
result, and a dashed line indicates the linear damping rate. The time is normalized by the 
inverse electron plasma frequency. 
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Figure 5: Electron phase space distributions (negative images with contours) at (a.b) 
t = 25 and (c,d) t = 50 in the two stream instability simulation, with (a,c) the MMA and 
(b,d) the CIP-CSL2 schemes. The numbers of grid points are 64 and 128 in the x and v 
directions. The time, position, and velocity are normalized by the inverse electron plasma 
frequency, Debye length, and thermal velocity, respectively. 
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Figure 6: Electron velocity space distribution at a; = 4.0 and i = 50 in the two stream 
instability simulation. A blue line is obtained from the MMA scheme with 64 grid points 
in the velocity space. Black solid and dashed lines are from the CIP-CSL2 scheme with 
128 and 64 points, respectively. The time, position, and velocity are normalized by the 
inverse electron plasma frequency, Debye length, and thermal velocity, respectively. 
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Figure 7: Temporal variation of (a) the total energy and (b) the entropy (normalized by 
their initial values) in the two stream instability simulation. Color lines are obtained from 
the MMA scheme with the different numbers of grid points in the velocity space (green = 
32, blue = 64, and red = 128), and black lines are from the CIP-CSL2 scheme (dashed = 
64 and solid = 128 points in the velocity space). The time is normalized by the inverse 
electron plasma frequency. 
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Figure 8: Relative error of the total energy at i = 75, as a function of the number of 
grid points in the velocity space. Triangles and diamonds are the results obtained from 
the MMA and CIP-CSL2 schemes. Symbols with the same colors are obtained from the 
simulations with the same memory usage. 
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Figure 9: Two-dimensional solid body rotation problem of a symmetric gaussian 
profile after (a,b,c) 50 and (d,e,f) 300 rotations, calculated with (a,d) the MMA, 
(b,e) the CIP-CSL2, and (c,f) the backsubstitution schemes. Contour levels are 
[0.01,0.03,0.1,0.25,0.5,0.7,0.9]. The numbers of grid points are 34 x 34 for the MMA, 
42 X 42 for the CIP-CSL2, and 84 x 84 for the backsubstitution. 
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Figure 10: Two-dimensional solid body rotation problem of an asymmetric gaussian profile. 
The format is same as Figure HI 
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Figure 11: Two-dimensional solid body rotation and advection problem of a symmetric 
gaussian profile. The format is same as Figure [9l 
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Figure 12: Temporal variation of the standard deviations (a) and (b) Oy in the two- 
dimensional rotation and advection problem shown in Figure 1111 Solid, dash, and dot- 
dashed lines are obtained from the MMA, CIP-CSL2, and backsubstitution schemes. 
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Figure 13: Errors of the two-dimensional solid body rotation and advection problem of a 
symmetric gaussian profile as a function of the grid size. Triangles, diamonds, and crosses 
are the results obtained from the MMA, CIP-CSL2, and backsubstitution schemes. A 
dashed line represents the third order accuracy. Symbols with the same colors are obtained 
from the simulations with the same memory usage. 
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Figure 14: Fourier spectra of perpendicular-propagating waves obtained from the electro- 
magnetic (a) Vlasov simulation and (b) PIC simulation until tUpet — 723.4. Horizontal 
and vertical axes are the wavcnumbcr and frequency normalized by the inverse electron 
gyro radius and the electron gyro frequency. Color contour shows the Fourier component 
of Ex (normalized by its maximum value) to the power of 0.15 (for illustration). From top 
to bottom, dashed lines represent the R-modc cutoff, upper hybrid, L-mode cutoff, and 
lower hybrid frequencies. Dot-dashed lines represent the dispersion relation of the light 
mode in vacuum. 
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Figure 15: Fourier spectrum of perpendicular-propagating waves obtained from the elec- 
tromagnetic Vlasov simulation until ujpet = 1447.4. Horizontal and vertical axes are the 
wavenumber and frequency normalized by the inverse proton gyro radius and the proton 
gyro frequency. Color contour shows the Fourier component of (normalized by its 
maximum value) to the power of 0.15 (for illustration). A dot-dashed line represents the 
dispersion relation of the Alfven wave. 
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Figure 16: One-dimensional Harris current sheet simulation with the MMA scheme, (a) 
The magnetic field (green lines) and density (black lines) distributions, (b) The pressure 
distribution. Green, black, and red lines are the magnetic, plasma, and total pressures, 
respectively. Solid lines are the simulation result and dashed lines are the initial. The grid 
size is A/j. 
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